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Abstract 

We use a physically-based crystal plasticity model to predict the yield strength of body-centered cubic (bcc) 
tungsten single crystals subjected to uniaxial loading. Our model captures the thermally-activated character 
of screw dislocation motion and full non-Schmid effects, both of which are known to play a critical role in 
bcc plasticity The model uses atomistic calculations as the sole source of constitutive information, with no 
parameter fitting of any kind to experimental data. Our results are in excellent agreement with experimental 
measurements of the yield stress as a function of temperature for a number of loading orientations. The 
validated methodology is then employed to calculate the temperature and strain-rate dependence of the 
yield strength for 231 crystallographic orientations within the standard stereographic triangle. We extract 
the strain-rate sensitivity of W crystals at different temperatures, and finish with the calculation of yield 
surfaces under biaxial loading conditions that can be used to define effective yield criteria for engineering 
design models. 

Keywords: Bcc crystal plasticity. Yield stress, Non-Schmid effects. Screw dislocations. Single crystal 
tungsten. Uniaxial/biaxial loading 


1. Background and motivation 

The plastic behavior of body-centered cubic (bcc) single crystals at low to medium homologous tempera¬ 
tures is governed by the motion of 11) screw dislocations on close-packed crystallographic planes. There 
are two particularities that make bcc metals unique in relation to their deformation characteristics. The first 
one is the thermally-activated nature of screw dislocation glide, a consequence of the compact (non-planar) 
structure of the dislocation core at the atomistic level (Vitek, 2004; Wurster et al., 2010; Li et al., 2012; 
Samolyuk et al., 2013). This feature is also responsible for the high intrinsic friction stresses reported in 
the literature for bcc metals and their alloys (Romaner et al., 2010; Samolyuk et al., 2013). The second is 
the breakdown of the standard geometric projection rule for the resolved shear stress (RSS) from the total 
stress tensor known as Schmid law (Schmid & Boas, 1935). This is owed to both specific crystallographic 
properties of the bcc lattice structure as well as to the coupling between the dislocation core and non-glide 
components of the stress tensor, which -to the best of our understanding- is unique to bcc crystals (Bulatov 
et al., 1999; Brinckmann et al., 2008; Woodward & Rao, 2001; Chaussidon et al., 2006; Groger & Vitek, 
2005). These anomalies have been the subject of much research and discussion going back to the 1960’s 
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(Takeuchi et al., 1967; Hull et al., 1967; Duesbery, 1969; Duesbery & Foxall, 1969), both experimentally 
and -more recently- using computational atomistic models. 

In regards to the first point above, at low stresses slip proceeds via the thermally activated nucleation of 
steps on the dislocation line, known as kink pairs, and their subsequent sideward relaxation. For a constant 
strain rate, this gives rise to the characteristic temperature dependence of the flow stress in bcc single 
crystals, which has been observed for all refractory metals and is considered to be a principal signature of 
their plastic response (Seeger, 1981; Ackermann et al., 1983; Taylor, 1992; Gordon et al., 2010; Chaussidon 
et al., 2006; Yang et al., 2001). The flow stress is considered to be composed of thermal and athermal 
contributions, with the latter depending on temperature only as the elastic moduli. Dislocation glide is 
thought to occur on {110}, {112}, and even {123} planes, depending on temperature and stress, over a 
periodic energy landscape known as the Peierls potential Up. The connection between the experimentally 
measured flow stress and this periodic energy potential is via the critical stress for which Up vanishes at zero 
temperature, known as the Peierls stress up. Theoretically then, the flow stress at very low temperatures 

25 K) is thought to represent the macroscopic equivalent of up as the temperature approaches 0 K. op 
can thus be unequivocally defined and has been the object of considerable numerical work since the first 
atomistic models were devised by Vitek and co-workers starting in the 1970s (Vitek & Yamaguchi, 1973). 

For their part, non-Schmid effects were detected in tests done in the 1930’s by Taylor in the wake of his 
seminal works on plastic flow and strain hardening (Taylor, 1928, 1934a,b). Subsequent observations and 
measurements (Sestak & Zarubova, 1965; Sherwood et al., 1967; Zwiesele & Diehl, 1979; Christian, 1983; 
Pichl, 2002; Escaig, 1968, 1974), and a rigorous theoretical formulation of the problem (Duesbery & Vitek, 
1998; Ito & Vitek, 2001; Woodward & Rao, 2001; Groger & Vitek, 2005; Chaussidon et al., 2006; Groger 
et al., 2008a,b; Soare, 2014) have established non-Schmid behavior as a principal tenet of bcc plasticity that 
must be accounted for in order to understand bcc plastic flow. In terms of phenomenology, the two essential 
aspects to bear in mind are (i) that the resolved shear stress is not independent of the sign of the stress 
in glide planes of the (111) zone (the so-called twinning/anti-twinning asymmetry), and (ii) that non-glide 
components of the stress tensor -i.e. those which are perpendicular to the Burgers vector- play a role on the 
magnitude and sign of the RSS on the glide plane of interest. 

Areas where we do not have a complete understanding of bcc plastic picture include the value of the 
flow stress at near zero absolute temperatures, the meaning of the so-called knee temperature, and the 
onset of athermal flow. In the last two decades, computer simulation has unquestionably emerged as 
discipline capable of shedding light on these processes on a similar footing with experiments, providing 
physically-substantiated explanations across a range of temporal and spatial scales. These include the use 
and application of density-functional theory methods (Ventelon & Willaime, 2007; Ventelon et al., 2013; 
Weinberger et al., 2013; Dezerald et al., 2014, 2015), semi empirical atomistic calculations and molecular 
dynamics calculations (Gilbert et al., 2011; Que 3 n:eau et al., 2011; Ghang et al., 2001; Komanduri et al., 
2001), kinetic Monte Garlo (Lin & Ghrzan, 1999; Gai et al., 2002; Deo & Srolovitz, 2002; Scarle et al., 2004; 
Stukowski et al., 2015), and crystal plasticity (GP)(Qin & Bassani, 1992; Dao & Asaro, 1993; Briinig, 1997), 
to name but a few. In general, while there is no doubt that the intricacies associated with ^(111) screw 
dislocation glide -including its thermally activated nature and deviations from Schmid law- cannot but be 
resolved using methods capable of atomistic resolution, one must recognize that, at the same time, flow is a 
phenomenon potentially involving statistically-significant amounts of dislocations and -as such- cannot be 
captured resorting to atomistic calculations only. 

Modeling thermally-activated flow and non-Schmid effects in bcc systems has been the subject of much 
work, starting in the 1980s and, particularly, in recent times. Different authors have considered different 
subsets of the {110}, {112}, and {123} families of glide planes, without (Raphanel & Van Houtte, 1985; 
Holscher et al., 1991, 1994; Raabe et al., 1994; Raabe, 1995a,b; Peeters et al., 2000; Stainier et al., 2002; 
Erieau & Rey, 2004; Ma et al., 2007; Hamelin et al., 2011; Kitayama et al., 2013) and with non-Schmid 
effects (Lee et al., 1999; Kuchnicki et al., 2008; Koester et al., 2012; Weinberger et al., 2012; Alankar et al., 
2014; Lim et al., 2013; Narayanan et al., 2014; Patra et al., 2014; Knezevic et al., 2014; Lim et al., 2015a,b). 
Of particular interest are some recent simulations where the flow rule is directly formulated on the basis 
of screw dislocation properties in Fe (Yalcinkaya et al., 2008; Koester et al., 2012; Alankar et al., 2014; 
Narayanan et al., 2014; Patra et al., 2014; Lim et al., 2015b), Ta (Kuchnicki et al., 2008; Lim et al., 2013; 
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Knezevic et al., 2014; Lim et al., 2015a), Mo (Yalcinkaya et al., 2008; Weinberger et al., 2012; Lim et al., 
2013, 2015a), W (Lee et al., 1999; Weinberger et al., 2012; Lim et al., 2013; Knezevic et aL, 2014; Lim et al., 
2015a), and Nb (Yalcinkaya et al., 2008; Lim et aL, 2015a). These works also include non-Schmid effects 
following the model proposed by Vitek and Bassani (Duesbery & Vitek, 1998; Qin & Bassani, 1992; Groger 
et al., 2008a,b). However, albeit very useful for certain applications, all these works resort to (i) a partial 
consideration of non-Schmid effects, and (ii) some kind or another of parameter fitting with experimental 
data, which prevents their use in regions of the parameter space outside the range of fitting and does not 
link the effective (macroscopic) response to exclusively fundamental material properties and features. 

In this work, we provide a unified computational methodology consisting of rate-dependent crystal plas¬ 
ticity calculations parameterized entirely and exclusively to atomistic calculations. We show that a full 
description of non-Schmid effects, together with the state of the art in terms of our understanding of 
thermally-activated screw dislocation motion, suffices to capture the experimentally measured tempera¬ 
ture dependence of the flow stress in tungsten. This is achieved in a fully classical framework, without the 
need for quantum effects recently invoked to explain the long standing discrepancy observed between the 
experimentally-measured flow stress below 25 K and calculated values of the Peierls stress (Proville et al., 
2012). Our methodology also captures the athermal limit of W to within 5% of the experimental value. We 
emphasize that this agreement is reached without fitting to any experimental data, all the parameterization 
is done from first principles atomistic calculations. 

Our paper is organized as follows. After this introduction, we provide an overview of the CP method in 
Section 2.2. This is followed by Sections 2.3.1 and 2.3.2, where the formulation of the dislocation mobility 
law and the implementation of non-Schmid effects are presented, including a detailed description of the 
parameterization procedure employed. The results are given in Section 3, which includes: (i) the validation 
exercise, with special focus on uniaxial tests as a function of temperature for several loading orientations; 
(ii) the calculation of temperature and strain rate dependence of the yield strength for uniaxial tensile 
tests as a function of orientation; and (iii) yield surfaces under biaxial loading conditions as a function of 
temperature. We finalize in Section 4 with a brief discussion and the conclusions. 


2. Computational methods 


2. 1 . Flow kinematics 

Roters et al. (2010) have presented a detailed review of the kinematic and constitutive aspects of crystal 
plasticity and here we simply provide a brief overview of the fundamental theory. The kinematics for elasto- 
plastic behavior is defined within the finite deformation framework. The material deformation involves both 
a reversible lattice response to externally imposed loads or displacements (elastic), and a permanent defor¬ 
mation (irreversible shape change) that remains after all external constraints cease to be applied (plastic). 
Consequently, crystal plasticity formulations rely on the definition of three reference systems: (i) a fixed 
coordinate system that represents a laboratory (undeformed) frame of reference, (ii) a current (also known 
as material) frame of reference that represents the global (deformed) shape of the material, and (iii) a lat¬ 
tice coordinate system that represents distortions of the underlying crystal structure of the deformed body. 
Although reference system (i) is used for mathematical convenience, the distinction between (ii) and (iii) 
is necessary to calculate internal stresses, which arise from distortions defined with respect to a crystallo¬ 
graphic frame of reference, as global shape changes may not necessarily have a one-to-one correspondence 
to internal lattice distortions (Lubliner, 2008; Roters et al., 2010). 

Mathematically, each point X in the reference configuration may be mapped to its image in the cur¬ 
rent configuration x by means of a linear transformation represented by the deformation gradient tensor F, 
defined as: 



( 1 ) 


In general, F is not a symmetric tensor. However, invariance requirementsmake it more desirable to work 
with symmetric measures of strain. One such measure is the so-called Lagrangian strain: 
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( 2 ) 


which refers the deformation of the solid to the reference configuration (J is the identity tensor). In the 
above equation, C is the so-called right Cauchy-Green tensor. 

Following Lee (1969), the total deformation gradient F can be multiplicatively decomposed into an 
elastic, F^, and a plastic, Fp, part^, i.e.\ 

F = FeFp (3) 

whence 

F,=FFy^Fp=F-^F 

This is schematically shown in Figure 1, where the relationship between the reference, intermediate, and 
current configurations is provided. To close the CP model, we take the time rate in eq. (1), which results in: 

L = FF-i = F,F-^ + F, (fpf-1) F-^ =L,+ F.L^FJ^ (4) 

where Lp is the plastic velocity gradient, which is evaluated in the intermediate configuration and must 
therefore be mapped into the current configuration by F^. Constitutive information enters the CP model via 
Lp, which is described in the following section. 



Figure 1: Multiplicative decomposition of the deformation gradient F. 

The above finite-deformation kinematic frameworkis implemented into the Diisseldorf Advanced Mate¬ 
rials Simulation Kit (DAMASK), which is the tool employed in this work to carry out of the calculations. 
DAMASK is a flexible and hierarchically structured model of material point behavior for the solution of 
elastoplastic boundary value problems along with damage and thermal physics (Roters et al., 2012). 

2.2. Solution procedure and constitutive model 

A Hookean constitutive response is assumed such that the stress depends linearly on the elastic strain via 
the anisotropic elastic stiffness tensor C. Both the stress and strain measures that are used internally are 
formulated in terms of material coordinates. For the stress, we use the second Piola-Kirchhojf stress measure 
S, defined as: 

(5) 


^ It must be noted that other decompositions are also admissible (Fish & Shek, 2000) . The reader is referred to the work by Reina & 
Conti (2014) for a discussion on the uniqueness and validity of the multiplicative decomposition. 
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where is the (elastic) Green-Lagrange strain tensor. S and are both symmetric material tensors, and 
thus C is itself symmetric such that a general 3 x 3 x 3 x 3 tensor can be written as a 6 x 6 matrix. For cubic 
lattices, C can be reduced by symmetry to only the three independent elastic constants Cn, C 12 , and C 44 . 

The stress S acts as the driving force for the plastic velocity gradient Lp. Lp depends on the underlying 
microstructure via a set of state variables ^ defined by the plasticity model employed: 

( 6 ) 

Lp controls the evolution of the plastic deformation gradient: 

Fp = LpEp (7) 

The set of nonlinear eqs. (3) and (5) to (7) must be solved iteratively, which in DAMASK is done by using 
an integration algorithm based on the implicit scheme originally proposed by Kalidindi et al. (1992). The 
linear system is solved iteratively using the Newton-Raphson technique and, once convergence is achieved, 
the plastic deformation gradient is obtained using the Euler backward update. In this integration scheme 
described above, the primary variable to solve for is the plastic velocity gradient. However, one may devise 
schemes where the primary variables are the stress, the plastic or elastic deformation gradient, the internal 
variables or a combination thereof. Such schemes may be chosen on the basis of computational efficiency 
(Dumoulin et al., 2009). 

Constitutive information for the plastic regime enters the CP model via eq. ( 6 ), where the dependencies 
of the flow rule on each of the state variables are established. It is here where the plastic deformation modes 
are defined, their geometric particularities, as well as specifics associated with the crystal structure under 
study. The CP model must also include evolution equations for the state variables $: 

( 8 ) 

where the details again depend on the model selected. In DAMASK, various integration schemes for the 
state update exist (Roters et al., 2012). Then, two integration schemes are performed staggered: eqs. (3) to 
(7) are solved at a fixed plastic state, followed by a state update. This procedure is iteratively repeated until 
a converged solution is achieved within the given tolerances. More details about the implementation of this 
technique in the code are given by Kalidindi et al. (1992). In general then, the stress in the CP model can be 
considered a response function of the position r, the deformation state E, the set of state variables and a 
set of boundary conditions, i.e. 

S=f{r,E,^,...) (9) 

In these calculations we are interested in simulating engineering stress-strain tests, and -consequently- 
it is helpful to express the results in the reference coordinate frame. For the stress, we use the first Piola- 
Kirchhoff measure defined as: 

P=FeS 

Note that, although in general P is not symmetric, for uniaxial tension tests P is a symmetric tensor on 
account of E^ being symmetric. For the strain, we use the Biot tensor: 

B = U-I = VF^-I 

where U is the right stretch tensor^. The stress-strain curves shown throughout this paper are obtained by 
tracking the evolution of P^z where z is designated as the loading direction. For uniaxial loading 

. 7 . . . 

simulations, E = E and thus C = E ^ B. We refer to the deformation rate represented by F^z genetically as 
£ in the remainder of this paper. 


emerges from the so-called polar decomposition: F = RU, where R is a matrix the represents rigid rotation, and U is a pure 
stretch. Plasticity-induced crystal rotations are very important and give rise to crystallographic texture evolution in deformed crystals. 
However, only yielding is of concern here, and thus U is the component of interest. 
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2.3. The flow rule 

In the present CP calculations it is assumed that all the plastic deformation is due to dislocation slip. 
Then, the plastic velocity gradient can be written as: 

h=l]Psr ( 10 ) 

a 

where is the slip rate on slip system a, and Pg is a geometric projection tensor that will be defined later. 
The slip rate is calculated from Orowan’s equation: 

y^ = bp^Vs{T^,T) ( 11 ) 

where b = uq a/3/2 is the modulus of the Burgers vector, uq is the lattice parameter, T the absolute temper¬ 
ature, is the (mobile) screw dislocation density in slip system a, and Vs{t^,T) is the screw dislocation 
velocity. The present formulation of the flow rule belongs to the class of non-associated, rate-dependent CP 
models (McDowell, 2008). 

The two characteristics that are particular to bcc plasticity are the thermally-activated nature of screw 
dislocation motion, which makes it the rate-controlling process during plastic deformation, and the exis¬ 
tence of non-Schmid effects, i.e. deviations from the geometric projection law for the resolved shear stress. 
Both of these physical processes have been known for several decades, and have been carefully analyzed 
experimentally (cf. Section 1). If our intent is to predict the temperature dependence of the flow stress in 
bcc metals, accurate physical descriptions of both of them must be incorporated into our CP model. This is 
the subject of the following sections. As we shall see, non-Schmid effects establish the form of the projection 
tensor Pf^^, while the velocity T) captures the thermally activated character of dislocation motion. We 
make tungsten the object of our study for a number reasons presented in previous works (Cereceda et al., 
2013; Stukowski et al., 2015)^. 


2.3.1. Screw dislocation mobility law 

Except at high homologous temperatures and strain rates, screw dislocation motion is the rate-limiting 
step in bcc crystal deformation. Although recent dislocation dynamics simulations in a-Fe challenge the 
notion that the dislocation density is monolithic across the entire temperature range (Monnet et al., 2009; 
Naamane et al., 2010; Monnet et al., 2011; Tang & Marian, 2014), it is reasonable to assume a dominance 
of screw dislocations in the temperature and strain rate regimes considered in this work (0 < T/T^ < 0.2 
and i ss 10“^ s“^). In the thermally activated regime, screw dislocation motion proceeds via the nucleation 
of kink-pairs and their subsequent lateral relaxation. In the regime of interest here, kink relaxation is a 
significantly faster process than kink-pair nucleation, and it can thus be assumed that no new kink-pairs will 
be nucleated while lateral kink motion is underway (Stukowski et al., 2015). Such assumption leads to the 
following expression for the total time, tf, required for a kink pair to form and sweep a rectilinear screw 
dislocation segment of length lying on a given slip plane: 


tf - tn + tk - /(t^, T) ^ + 


-w 

2vk{r^,T) 


( 12 ) 


where is the mean time to nucleate a kink pair, tk is the time needed for a kink to sweep half a segment 
length, 7 is the kink-pair nucleation rate, w is the kink-pair separation, and Vk is the kink velocity. The 
kink-pair nucleation rate follows an Arrhenius formulation (Stukowski et al., 2015): 


/(t^T) 


b 



r J 


( 13 ) 


is an elastic isotropic metal, which simplifies the constitutive plastic formulation. 
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where Vq is an attempt frequency, AH^p is the activation enthalpy of a kink pair at stress and k is 
Boltzmann’s constant. For its part, the kink velocity can be expressed as (Dorn & Rajnak, 1964; Kocks et ak, 
1975): 

nir^T)^^) (14) 

where B is friction coefficient typically assumed to be linearly dependent on temperature. However, calcu¬ 
lations made to obtain the value of B for the interatomic potential employed in this work, have yielded no 
temperature dependence, and here B is simply a constant (Swinburne et al., 2013) . The dislocation velocity 
can be obtained after operating with eqs. (13) and (14) as: 


h h 
h tn + h 


2bhx^VQ{X^ — w) exp 

2b^x^ + vqB(A^ — ipyexp 


(15) 


where h = uq a/6/3 is the distance between two consecutive Peierls valleys. We note that at low temperatures, 
or when tj^ « t^, the second term in the denominator vanishes and one recovers the standard diffusive 
velocity expression commonly used in crystal plasticity and dislocation dynamics: 





w) 


exp 





Figure 2: Schematic depiction of a kink pair on a screw segment of length A lying on a slip plane (of the {110} family). The vertical 
axis represents the potential energy with the Peierls potential clearly marked. The dashed line represents the initial equilibrium line 
position. 

The parameterization of eq. (15) is a critical step that establishes a physical connection with the scales 
where kink-pairs are resolved as atomistic entities. This is the first essential piece of physics required to 
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achieve predictive capabilities. We have devoted much effort in past works to calculate the necessary pa¬ 
rameters from fundamental models based on semiempirical interatomic potentials (Cereceda et ah, 2013; 
Stukowski et ah, 2015). The list of parameters employed in this work and their associated values and units 
are given in Table 1. The physical meaning of some of these parameters is best expressed in pictorial form. 
Figure 2 shows a schematic diagram of the topology of a kink pair lying on the Peierls energy substrate. The 
figure highlights the physical meaning of each parameter listed in the table. In addition to the references 
provided earlier, a detailed description of the protocols used to calculate all the adjustable parameters in our 
formulation is provided by Cereceda (2015). 

At this stage, it is worth to introduce a note about the available slip systems (which establish the running 
indices of a, Stukowski et al. (2015) have shown that in W an elementary glide on a {112} plane is 
a composite of two elementary steps on alternate {110} planes. Judging by these results, we conclude 
that glide on any given plane is achieved by way of sequential { 110 } jumps, which constitutes the basis 
to simulate plastic yielding in the foregoing Sections. This is consistent with recent atomistic simulations 
(Cereceda et al., 2013) and experiments (Caillard, 2010a,b; Marichal et al., 2013, 2014) and limits the 
number of available slip systems in our study to 12 (listed in Appendix A). We note that this model of slip 
for W is not necessarily suggestive of what may happen in other bcc crystals (Franciosi et al., 2015). 

2.3.2. Projection tensor and non-Schmid ejfects 

The tensor Pg introduced in eq. (10) represents the Schmid (geometric) projection of the strain rate 
contribution from a slip system defined by the plane normal and slip direction (both unit vectors). 
However, as pointed out above, Pg does not capture the full panoply of non-Schmid effects needed to 
calculate the value of the resolved shear stress on that slip system, r^. For this, we introduce a total 
projection tensor Pf^^ such that: 


_ TiOC 

~ ^ tot 


:a=(p«+P“/AT + P“g):a 


(16) 


where 

is the Schmid tensor, with 
the Cauchy (true) stress and / 


Pg =m^ 

a = J-^FSF'^ 

det(P) the Jacobian. The tensors 


Kg =a2{n^ xm^)(g)n^ + xm^)(x)n? 


(17) 


(18) 

(19) 


are non-Schmid tensors representing respectively the twinning/anti-twinning asymmetry (T/AT) and the 
effects due to non-glide stress components, ai, a 2 , and ^3 are material-dependent constants that must also 
be calculated and added to our parameterization database. The vector forms an angle of —60° with the 
reference slip plane defined by and changes sign with the direction of slip on each glide plane (Koester 
et al., 2012 ). 

The present non-Schmid formulation was originally developed by Vitek and expanded by others, and has 
been successfully used to propose yielding criteria adapted to finite element and crystal plasticity calculations 
in a number of cases (Groger & Vitek, 2008; Chen et al., 2013; Weinberger et al., 2012). The reader is 
referred to these works for more details but it is worth pointing out that the methodology that these authors 
have proposed is not unique, and that other rigorous implementations of non-Schmid effects could equally 
be devised. For the purposes of this section, suffice it to say that the particularities of the screw dislocation 
core and the bcc lattice structure result in deviations from a purely geometric projection. These deviations 
originate, respectively, from a geometric asymmetry between the twinning and anti-twinning directions of 
the ( 111 ) zone -from which ui is first calculated-, and from the effect that nonglide components (termed 
generically ‘a’) of the local stress tensor have on the critical resolved shear stress, from which ^2 
are obtained. Atomistic calculations specifically designed to calculate the non-Schmid critical stress 
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as a function of the angle x between the maximum resolved shear stress (MRSS) plane were performed 
according to the geometry shown schematically in Figure 3. The Figure shows the mapping between the 
atomistic box and the crystallography of the [111] zone. Following the sign convention used in the Figure, 
the stress tensor applied is: 

/ -d 0 0 \ 

Out (20) 

\ 0 T 0 / 

which activates axial (nonglide) stress components while maintainign zero pressure, is expressed as a 



Figure 3: Crystallographic diagram of the [111] zone in the bcc lattice with each {110} and {112} clearly labeled. The picture also 
shows a mapping of the [111] zone to a schematic atomistic box containing a screw dislocation subjected to shear and nonglide stresses 
according to Vitek’s convention. This setup is used to calculate the critical RSS using atomistic calculations (cf. Appendix B). The glide 
n^, auxiliary rii^ and MRSS planes are labeled in each case. A [101] glide plane corresponds to a = 2 in our CP calculations. 

combination of the contributions displayed in Fig. 3: 

'r;+a(a2sin(2;t)+a3sin(2;i: + f)) 

Tc — / Tr\ 

COSX + Cll COS [X + f j 

where r* is a fitting constant that represents the Peierls stress. The details of these atomistic calculations 
are provided in Appendix B. The results for are shown in Figure 4 as a function of x and a, with t*, ai, 
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a 2 , and ^3 given in Table 1. It is worth noting that the relation between and a has been established for 
tensile nonglide stresses only (n > 0), for consistency with the linear dependence used in the work of Vitek 
and collaborators (Groger et ah, 2008a, b) that has been used in other crystal plasticity works (Koester et ah, 
2012). However, nothing precludes the use of nonlinear fitting functions that capture both the tensile and 
compressive regimes simultaneously (cf. Appendix B). It is worth noting that Groger et al. (2008b) obtained 
values of ai = 0, a 2 = 0.56, and ^3 = 0.75 using a bond-order potential, substantially far from our values for 
those parameters. 



-30 -20 -10 0 10 20 30 

Z[°] 

Figure 4: Critical resolved shear stress as a function of the angle x between the MRSS and glide planes and the value of the nonglide 
stress component a with the sign convention according to Fig. 3. The value of the Peierls stress ap = 2.03 GPa is circled. 

By way of example, we calculate the maximum projection factor M for directions in the standard stere¬ 
ographic triangle using the fully parameterized projection tensor: 

M = {1®1): Ptot = {PI + P?/at + ^ng) (22) 

where I is the loading direction, which is obtained by visiting each of the nodes resulting from the discretiza¬ 
tion of the standard triangle area into a uniform grid consisting of 231 points. The results for tension (a > 0) 
are shown in Figure 5. It is clear than non-Schmid effects -particularly the impact of nonglide components- 
are critical to calculate the RSS on a given slip system. We find that from a maximum nominal value of 
M = 0.5 for the standard Schmid law there is a twofold amplification when the twinning/anti¬ 
twinning asymmetry is considered (P^^^ + and an astonishing fourfold increase when nonglide 

effects are also included (P^^^ + As we shall see in Section 3.2, this has extraordinary 

importance when comparing GP calculations to experimental measurements. 
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omax I omax i omax 
+ ^T/AT + ^ng 



[10 1 ] 


Figure 5: Projection factor according to eq. (22) for 231 directions within the standard triangle. The contributions of each of therms in 
eq. (16) are broken down for comparison. 


2.4. Dislocation density evolution model 

To close the model, one needs to provide an evolution law for the dislocation density in Orowan’s equa¬ 
tion 11. There are numerous density evolution models proposed in the literature, each with a specific 
domain of applicability (Mecking & Kocks, 1981; Estrin, 1996; Arsenlis & Parks, 2002; Stainier et ah, 2002; 
Barton et ah, 2011). In this work we are mainly interested in yielding, i.e. the elastic-to-plastic transition 
before dislocation-based slip takes on a dominant role in the constitutive model. We use the model presented 
by Roters (2011), in which the mobile dislocation density on slip system a evolves in time according to: 


P“=PLlt + Pann (23) 

The evolution model is initialized by the dislocation density at t = 0, pQ. In eq. (23), p^^^^ and p^^n represent 
the dislocation multiplication and dislocation annihilation rate terms, respectively. In this model, both p^^^^ 
and p^nn directly proportional to the plastic strain rate. Dislocation multiplication is treated as being 
proportional to the inverse mean free path of the dislocations, : 


■a 

rmult 


(24) 


11 





which is defined as a function of the grain size dg, the forest dislocation density pj, and a hardening constant 


c: 


j_ _ j_ 

A“ dg 


Here, c and dg are set, respectively, to one and to an arbitrarily high value such that the term controlling the 
dislocation mean free path is: 

The forest dislocation density is calculated as (Roters et ah, 2010): 


pJ=Y,P^\n- 


m' 


(26) 


Note that, in general, the mean free path as defined in eq. (25) need not be equal to the effective dislo¬ 
cation segment length (in fact, it can be up to several orders of magnitude different (Basinski & Basinski, 
1979)). However, our model is designed with well-annealed, high-purity single W crystals in mind, with 
low initial dislocation densities and no impurities or obstacles other than dislocations themselves. Under 
this assumption, the use eq. (25) can be justified in this case (Stainier et al., 2002; Monnet et al., 2011). 

For its part, dislocation annihilation occurs spontaneously when dipoles approach to within a spacing of 


d 


edge* 




(27) 


Equations (23) through (27) form the basis of the Kocks-Mecking family of dislocation density evolution 
models (Mecking & Kocks, 1981). These models offer two interesting connections with the broader CP 
formulation employed here. First, a relation between the dislocation density evolution model and Section 
2.3.1 is established by way of the dislocation mean free path A^, which determines the available segment 
length in the dislocation mobility function (eq. (15)). In this fashion, the dislocation velocity -and, through 
it, the plastic strain rate- is self-consistently linked to the microstructure changes predicted by the model. 
Second, by virtue of the existence of latent and self-hardening, the model provides a correction to the 
available RSS for dislocation motion in eq. (15) of the following form: 


- t,, = P“ t :cT-pb 

where is the hardening stress and are the coefficients of the interaction matrix, which character¬ 
izes the interaction strength between slip systems a and n’ as a result of six possible independent interac¬ 
tions (Franciosi, 1983, 1985): self, coplanar, collinear, mixed-asymmetrical junction (orthogonal), mixed- 
symmetrical junction (glissile) and edge junction (sessile) (Madec & Kubin, 2004). The values of 
employed here are given in Table 2, and were obtained from dislocation dynamics simulations of isotropic 
elastic bcc Fe under uniaxial deformation^ (Queyreau et al., 2009). The correspondence between each co¬ 
efficient and each slip system considered in this work is given in Appendix A. replaces in eqs. (12) 
to (15), although, as mentioned earlier, this pertains mainly to the plastic flow regime and -as such- is not 
expected to have a significant bearing on our calculations of Uy. 


3. Results 

In this Section we present results of uniaxial and biaxial tensile test simulations to explore the depen¬ 
dence of the yield strength on loading direction, temperature and strain rate. First, however, a robust and 
consistent yield criterion must be defined to extract the yield stress from the raw output data from DAMASK. 

^Although the coefficients were calculated for bcc Fe and not W, the results are equally applicable because Fe was treated 
as isotropic elastic -as is W- and the interaction matrix coefficients are non-dimensional and independent of the value of the plastic 
constants considered. 
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Table 1: List of parameters and functional dependences for fitting the CP model. All of these parameters have been obtained using 
dedicated atomistic calculations. The parameter s represents the normalized shear stress: s = ^ (cf. eq. (28)). 


parameter 

value or function 

units 

Uq 

3.143 

0 

A 

b 

2.72 

A 

h 

flO V^/S 

0 

A 

Cii 

523 

GPa 

Ci2 

202 

GPa 

C 44 

161 

GPa 

Vo 

9.1 X IQii 

s“^ 

ap 

2.03 

GPa 

B 

8.3 X 10-5 

Pa-s 

AH{s;T) 

AHo(l-sP)‘? 

eV 

AHo 

1.63 

eV 

P 

0.86 

- 


1.69 

- 

w 

11 

b 


T* +cr(fl2sin(2x)+fl3 sin(2x+7t/6)) 

GPa 


cosx+ai cos(7t/3+x) 

Cl\ 

0.938 

- 

a 2 

0.71 

- 


4.43 

- 

T* 

2.92 

GPa 

C 

1 

- 

dg 

2.72 

0 

A 

^edge 

2.72 

A 


Table 2: Values of for latent hardening in bcc crystals (from Queyreau et al. (2009)). 


self 

coplanar 

collinear 

orthogonal 

glissile 

sessile 

0.009 

0.009 

0.72 

0.05 

0.09 

0.06 


3.1. Yield criterion 

In metals, where dislocation flow is not a singular event but a diffuse continuous process, it is gener¬ 
ally accepted that the definition of yield point^ is not unique. Perhaps as the result of these conceptual 
indetermination, modern usage has evolved into that of an arbitrary rule, the 0 . 2 % strain offset rule for 
obtaining the yield stress of metals. For materials having nonlinear elastic behavior, there are not even 
arbitrary rules, only individual preferences and proclivities in defining yield when a given amount of strain 
has been reached. It is quite apparent then that to define robust yield criteria it is necessary that they be 
implemented and supported by consistent and meaningful definitions in terms of the stress-strain behavior. 
This is often difficult when the transition from the elastic to the inelastic regimes is obscured in the global 
picture of deformation. However, in the present calculations we effectively possess an arbitrary degree of 
data resolution and can define an unambiguous mathematical criterion. 

The preferred method for defining the elastic limit of a ductile material is to compute the second deriva¬ 
tive of the stress-strain curve, referred to genetically as (j{e), and identify the location of the inflection 


^Also referred to as elastic limit, proportionality limit, yield stress, etc. 
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point (Christensen, 2008). The yield point then corresponds to the strain, Cy, for which 
Mathematically: 

Gy = G{£y), ey:=e\ max 


d^G 




is maximum. 


(29) 


For ductile metals, the location of the maximum of the second derivative represents the point at which 
dislocation-mediated flow is the major contribution to L (cf. Section 2.1). However, this condition works 
surprisingly well for other materials such as glassy polymers, where flow might be caused by molecular 
rearrangement and damage at both the molecular and macroscopic scales (Bowden & Jukes, 1972). 

To illustrate the accuracy of the second-derivative method, we plot in Figure 6 the first and second 
derivative of a stress-strain curve corresponding to a [101] uniaxial tensile test of a W single crystal under 
representative initial conditions. Recall from Section 2.2 that the stress and strain metrics of choice are P 
and B, and so we plot 4^ and specifically. The inflection point -marked by a vertical dashed line in 

dBzz 



0.00 0.10 0.20 0.30 0.40 

B,, [%] 


Figure 6: Evolution of the stress Pzz with deformation B^z during a CP simulation of a uniaxial tensile test with [101] loading 
orientation (as depicted in the standard triangle). The first and second derivatives of the stress wr.t. to the strain are also plotted to 
illustrate the method of identification of the yield point according to this criterion. Also shown is the intercept of the curve with the 
0.2% strain offset criterion line. 

the figure- occurs for £y = 0.1105%, for which a value of Gy = 0.452 GPa is obtained. The figure also shows 
the 0.2% strain offset criterion, which -by contrast- gives £y = 0.3167% and Gy = 0.479 GPa, i.e. a three-fold 
difference in strain and approximately a 6% difference in stress with respect to the stress second derivative 
criterion. 

However, determining the first and second derivatives of the stress-strain relation can become numeri¬ 
cally intensive, especially when evaluating thousands of curves as is the case in this work. An approximation 
to this method that works particularly well for linear-elastic materials that display a clear elastic-to-plastic 
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transition is to take the yield point as the first point in the a(£) function that satisfies: 


da 

de 


<£( 1 - 6 ) 


i.e. Gy is measured as the stress for which a departure from linearity (as set by the elastic regime) larger 
than some small value 6 is observed in the stress-strain relation. We have found that a value of 6 % 0.01 
is sufficient to predict the value of Oy within a small error relative to the value furnished by the second- 
derivative method. By way of example, for the curve shown in Fig. 6 and 6 = 0.01, we find a values of 
Ey = 0.1055% and Oy = 0.435 GPa, or less than a 4% difference with the numbers according to the second- 
derivative criterion. With this reasonable accuracy and the computational advantages alluded to above, we 
then use the 6 = 0.01 criterion in the remainder of this paper. 


3.2. Model validation and initial results 

Prior to deploying our fully-parameterized CP method for numerically-intensive calculations, it is essen¬ 
tial to undergo a thorough exercise of validation. Experimental data from uniaxial tensile tests in single 
crystal W at low strain rates are scant and sporadic, with the main sources listed below: 

1. Argon & Maloof (1966) performed some early experiments at a strain rate of 10“^ s“^ and tempera¬ 
tures of 77, 199, 293, 373, and 450 K. These authors measured the yield strength for the three vertices 
of the stereographic triangle [001], [110], and [111] with an initial dislocation density of po ^ 10^^ 
m~^. 

2. Raffo (1969) analyzed the yielding behavior of arc-melted W between 77 and 680 K at £ = 8.3 x 10“^ 
s“^. However, the loading orientation is not given and most of the tests were done in compression. 

3. Stephens (1970) has carried out compression tests at 150, 300, and 590 K. This researcher focuses on 
dislocation density evolution and dislocation substructures, however, with a value ofpo ^ 1.4x 10^^ 
m“^, notably larger than in other tests. There have been other works that have also focused mainly 
on compression tests (Klopp et al., 1964; Gupta & Li, 1970). 

4. Brunner (2000, 2010) has performed a series of experiments more recently at temperatures between 
77 and 800 K. They employed a value of £ = 8.5 x 10“^ s“^ and loaded the system uniaxially along 
the [14 9] direction with a starting dislocation density of 5.5 x 10^ m~^. 

As pointed out in Section 2.3.2, our CP model is parameterized for tensile tests only and so for validation 
we focus on the works by Argon & Maloof (1966) and Brunner (2000, 2010). Argon & Maloof (1966) 
centered on multislip by considering mainly loading orientations coincident with the vertices of the standard 
triangle. Consequently, we replicate their test conditions in our CP model and compare the results obtained 
by taking into account all the different elements of the projection tensor (16). The results are shown in 
Figure 7 for the [111] and the [110] loading orientations, with the insets in both figures showing the relative 
importance of considering each of the non-Schmid contribution to the projection tensor incrementally. While 
our calculations are in general good agreement with the [111] test data, they deviate from the experimental 
results at the two lower temperature points for the [110] orientation. Argon & Maloof (1966) point out that, 
at low temperatures, deformation by twinning may play a larger role when loading along [110] relative to 
other orientations. This may be at the origin of the discrepancy, as twinning is not part of the catalog of 
deformation mechanisms considered in this model. 

Next we simulate uniaxial tensile tests under single slip conditions, i.e. along crystal orientations near 
the center of the standard triangle. This corresponds to the experiments by Brunner (2000, 2010) referred 
to above, which were done more recently with more advanced instrumentation. The results are shown in 
Figure 8, where we also show the curves using the different elements of eq. (16). This time, the agreement 
is striking, particularly again at temperatures above 400 K. Specifically, the athermal limit (% 710 K) is 
particularly well reproduced, as is the extrapolated critical stress at 0 K (Peierls stress), which is within 10% 
of the experimental values. 

Although, as noted earlier, the main focus of this work is on yielding, we have applied the fully param¬ 
eterized model to study the flow stress regime for some selected cases in Appendix C. The results shown 
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there demonstrate the performance of the method outside the primary range of application. While the 
model cannot be assumed to be predictive in the post-yield regime under general loading conditions, these 
are encouraging results that strengthen the notion that parameter-free CP calculations can perform well 
under specific deformation scenarios. 

With the confidence conferred on our CP model by the validation procedure, next we proceed to calcu¬ 
late the yield strength for a number of numerically-intensive scenarios. This is the object of the following 
sections. 

3.3. Uniaxial tensile tests 

In this Section, we report on the uniaxial yielding results as a function of temperature and strain rate. 
Our results are organized by strain rate, such that we first provide a detailed account of all the calculations 
at a given strain rate followed by a study on the dependence with L 

3.3.1. Results at e = 10“^ s“^ 

For these calculations, we have discretized the area of the standard triangle into a uniform grid consisting 
of 231 nodes, each representing a crystallographic loading orientation. We begin with calculations at a 
prescribed strain rate of £ = 10“^ s“^. Figure 9 shows colored contour plots of the yield stress in the 
lOO-to-600 K temperature range. Areas of high relative yield strength can be seen to concentrate around 
the vertices of the standard triangle, representing multislip conditions, whereas soft regions develop in two 
distinct locations of the triangle, one near the [324] zonal axis that then rotates towards [112] above 500 
K, and another near [102]. Note that, to accentuate the differences between hard and soft regions, each 
contour plot has its own specific numerical scale. 

We have extracted the specific location of the global extrema in the standard triangle and plot it as a 
function of temperature in Figure 10(a). The hardest direction is consistently the [101], while the softest is 
seen to revolve around the vicinity of the [112] axis, first along [30 18 41] at 100 K, then along [180 131 271] 
between 200 and 500 K, and finally rotating towards [9 9 34] for T > 500 K. Next, we plot the detailed 
temperature dependence of the yield strength corresponding to the hardest and softest directions -as given 
by Fig. 10(a)- for this strain rate in Figure 11. As the calculated data show, there is approximately a 30% 
difference in yield stress between the hardest and softest directions. Interestingly, this gap appears to be 
fairly independent of temperature. Above 650 K, the curves begin to level off, signaling the onset of the 
athermal regime. 

3.3.2. Dependence on strain rate and strain rate sensitivity 

In this Section we expand the analysis presented in the previous Section to strain rates of 10“^ and 
10“^ s“^. To avoid redundancies, here we show only the temperature trajectory of the softest and hardest 
loading orientations in Figs. 10(b) and 10(c), which emanate from calculations as those presented in Fig. 
9. The results are quantitative similar to the case of £ = 10“^ s“\ with the only appreciable deviations 
occurring at temperatures above 450 K. At these high temperatures, the softest orientation rotates clearly 
towards the vicinity of the [113] zonal axis, without excursions near [103] as was the case for the i: = 10“^ 
calculations. 

As above, we add the temperature dependence of the yield stress for the hardest and softest directions 
at these strain rates to Figure 11. The data show the same qualitative trend for all strain rates, with the 
same approximate 30% difference between the hard and soft orientations. However, useful information 
can be extracted if the strain-rate dependence of the yield stress is plotted for selected orientations. Then, 
one can calculate the so-called strain rate sensitivity, characterized by the strain rate sensitivity exponent m, 
of the material as a function of temperature. Strain rate sensitivity is exceedingly important to delay the 
onset of inhomogeneous deformation (Hutchinson & Neale, 1977), e.g. necking, and is used as a criterion 
to assess the possibility of superplastic behavior in certain kinds of materials (Hedworth & Stowell, 1971; 
Ariel! & Rosen, 1976). This belongs more in the realm of failure and is thus outside the scope of this paper. 
However, it is of interest to calculate the strain rate sensitivity of the yield stress and relate our findings to 
the larger failure picture if possible. 
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This precisely what is done in Figure 12 for [101] loading tests. The figure shows the variation of the 
yield strength at the three strain rates considered here, again in the range 100 < T < 600 K. The data can 
then be fitted to the following expression: 

(jy = (30) 

where C is a fitting constant. The strain rate sensitivity exponent is formally defined as: 


aiogUy 

dloge 


(31) 


m is plotted in the inset to Fig. 12, where it can be seen that it increases monotonically with temperature 
from a value of m = 0.01 at 100 K to % 0.2 at 600 K. The implications of these results will be discussed in 
Section 4. 


3.4. Biaxial loading tests and yield surfaces 

For non-associated CP formulations such as the present one, yielding is not a separate and independent 
criterion, but a consequence of the constitutive law of the material behavior (Bodner, 1968). Indeed, with 
yielding defined on the basis of the identification criterion introduced in Section 3.1, yield surfaces are 
furnished as a product of the CP calculations. In this Section we calculate the yield curves under biaxial 
stress conditions for selected pairs of orthogonal loading directions ly and As noted in Section 2.3.2, the 
present implementation of the non-Schmid stress projection law is only valid for tensile conditions^. Thus, 
our yield curves are only meaningful in the positive stress quadrant (or octant, for yield surfaces). The 
procedure to calculate each point of the yield surface consists of deforming the system simultaneously along 
the prescribed orientations until the material yields on either one according to criterion (29). The stresses 
Pzz and Pyy are then measured along both directions and the resulting duplet is added to the curve. Plane 
stress conditions are adopted along the remaining direction, i.e. P^^ = 0. The calculations are done at a 
nominal strain rate of e = s“\ with slight variations above and below this value in one of the loading 

directions to accumulate different levels of stress and map the entire stress quadrant. 

First we calculate the yield curve for ly = [111] and =[ll 2] as a function of temperature. Results are 
shown in Figure 13. The curves enclose domains that are everywhere convex, thus satisfying the Drucker- 
Prager criterion for stable plastic flow materials (Prager, 1952; Drucker et al., 1952). The absolute values 
and the temperature sensitivity of the yield stresses for the end cases of P^z = 0 Pyy = 0 are consistent 
with the results shown in Section 3.3 for the ly and chosen here. 

The next series of calculations involves determining the entire yield surface of the [111] zone, i.e. for a 
set of directions orthogonal to [111] in 10° intervals, at a fixed temperature of 300 K. Results are shown in 
Figure 14. Symmetry considerations limit the angular range to be explored to a 60° arc, which is shown 
in the figure. Yield surfaces such as this one are the culmination of crystal plasticity calculations, and 
can be used as constitutive input into continuum models to simulate effective mechanical behavior at the 
engineering scale, for component design and/or to simulate, e.g., thermo-mechanical treatments (Sheng 
et al., 2004; Serenelli et al., 2010). 


4. Discussion and conclusions 

In this Section we consider the most important implications of our results. First, we discuss one of 
the most salient characteristic of the current work. The present CP model uses a standard rate-dependent, 
finite-deformation, non-associated theory of crystal plasticity. However, while the underlying kinematic 
formulation serves as the mathematical framework upon which to build a physical methodology, it is via the 
connection to the material physics that the model is rendered truly predictive. Our technique does so by 
incorporating the following three features of bcc slip: 


^Although this is not a limitation in a strict sense as it is done simply for consistency with non-Schmid treatments published in the 
literature. 
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• A complete (T/AT plus nonglide) treatment of non-Schmid effects. 

• A kinematic flow rule based on a thermally-activated screw dislocation mobility. 

• Using accurate interatomic potentials for computing all the free parameters in the model. 

We have shown that the full model is capable of predicting the experimentally-measured temperature de¬ 
pendence of yield strength in the entire temperature range for W single crystals without parameter-fitting 
of any kind^. The sole source of material (constitutive) information is a carefully selected semi-empirical 
interatomic potential fitted exclusively to a DFT-generated dataset that includes the Peierls stress in its full 
atomistic meaning. This closes the gap seemingly separating electronic structure calculations of fundamental 
dislocation core properties and real measurements of the yield stress in uniaxial tensile tests of bcc materials. 

Indeed, much effort has been devoted to the study of this long-standing experiment/simulation discrep¬ 
ancy particularly at temperatures < 20 K. Explanations based on collective dislocation dynamics, such as 
network kinetics (Bulatov & Cai, 2002) and/or mutually interacting dislocations (Groger & Vitek, 2007) 
can be more or less discounted in light of recent detailed electron microscopy experiments of isolated screw 
dislocation motion (Caillard, 2010a,b, 2014). A more recent description, based on quantum effects at very 
low temperatures, has been put forward with reasonable success (Proville et al., 2012). On this basis, our 
first partial conclusion is that, while the present calculations do not provide sufficient grounds to invalidate 
these theories, they do clearly demonstrate that models based solely on classical mechanics -and without re¬ 
course to fitting to experimental results- can be formulated to predict the temperature dependence of the yield 
strength of bcc single crystals. Evidently, we issue this conclusion with caution, as W does not constitute by 
itself a representative sample to convincingly claim generality, but we believe that it constitutes a step in 
that direction. 

Another important physical aspect of tensile deformation in single bcc crystals is the seemingly distinct 
slip mechanisms operating in different temperature ranges. According to Seeger and collaborators, there 
are three clearly distinguishable temperature regions in the flow stress-temperature curves for bcc metals 
(Seeger, 1981, 1995; Brunner, 2000), namely, the so-called upper and lower bend temperatures, f and f, 
and the knee temperature t, f, and Tj^ delimit three different regimes where slip may occur on {110}, 
as well as {112}, glide planes, and give rise to different deformation mechanisms. Although these theories 
are substantiated by ample experimental data, there are recent studies that indicate that {110} slip may be 
sufficient to explain the most salient features of bcc plasticity (Marichal et al., 2013; Ali et al., 2011). This is 
consistent with the analysis presented here, backed by atomistic input, which suggests that only {110} slip 
is admissible in bcc W. Interestingly, the screw dislocation mobility law employed in this work, where {112} 
slip is disallowed by construction (cf. Section 2.3.1), is sufficient to quantitatively characterize the evolution 
of the yield stress across the entire temperature spectrum, without any ad hoc partition of mechanisms 
into different temperature regimes. We emphasize once more that the screw mobility law has been fitted 
exclusively to first-principles data. 

In Section 3.3.2 we have provided calculations of the strain rate sensitivity defined asm = dlogay/dlogL 
It must be noted that our value of m = 0.023 at 300 K obtained in the 10“^ > £ > 10“^ s“^ range is 
consistent with measurements performed by Zurek & Gray III (1991) in W compressed uniaxially at strain 
rates from 10“^ to 10^ s“^. Notwithstanding the differences in experimental methodology and strain rate 
regime, this is also encouraging agreement for a result other than yield, m is an important parameter 
for calculating the kink-pair activation enthalpy and activation volume from stress-relaxation tests. Note 
that some authors use an alternative definition for the strain rate sensitivity (Raffo, 1969; Brunner, 2000), 
namely, A = da/dloge, which is related to m via A = ma. We can then conclude that the agreement achieved 


^Of course, interatomic potentials -which form the basis of the constitutive information employed here- are subjected to a fair 
amount of fitting themselves, both to experimental data and first-principles calculations. However, potential fitting is extraneous to our 
work, in the sense that it was neither performed by us nor done with this application in mind, while the parameter fitting that we refer 
to here is dedicated specifically to reproduce experimental data of interest to the application of the model. 

^ T]^ is understood as the temperature above which the contribution of the kink-pair formation mechanism to the flow stress becomes 
negligibly small, i.e. it signals the athermal limit. 
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for a derivative quantity of the yield stress such as m is symptomatic of the quality of the method outside the 
primary validation space. 

The advantages of this and other CP methodologies w.r.t. more accurate techniques such as molecular dy¬ 
namics, dislocation dynamics, or phase field methods is of course their computational expediency Backed by 
the encouraging outcome of the validation exercise, this has enabled us to map the entire loading orientation 
space in the standard triangle (231 directions) as a function of temperature in a experimentally-meaningful 
strain rate range. These results can then be used to extract useful information, such as the strongest and 
softest orientations as a function of temperature and strain rate, or the strain rate sensitivity of our W model 
system. This information can ultimately be used to define yield criteria under a variety of conditions for 
more homogenized methods, with the aim put on component design. 

In this sense, the culmination of the CP simulations is the calculation of yield curves and yield surfaces 
in stress space. The stress space that we have chosen for our yield surface calculations is a purely biaxial 
one (in plane stress) with one fixed direction, chosen arbitrarily to be [111], and the family of orthogonal 
directions taken in 10° intervals. This biaxial loading configuration is the elementary basis for pressurized 
cylinders, e.g. pipes, and is thus useful to design components based on this geometry. As well, it can serve 
as the design premise for loaded plates under plane stress conditions. It is of interest to note that yield 
surfaces can also serve as the plastic potential in the fundamental theory of plasticity (Lubliner, 2008). This 
equivalence is valid when the critical resolved shear stress is not dependent on the current stress state^ 
(Lubliner, 2008; Starovoitov & Naghiyev, 2012). However, this may not be applicable in the present model, 
where the CRSS is seen to display a strong dependence on hydrostatic (nonglide) stress components as 
discussed in Section 2.3.2. This is also the case in rock and soil plasticity (e.g. Pariseau et al. (1968)). In 
such cases, the normality rule is referred to the pressure-dependent yield surface instead. 

A standing limitation of our model is that we have only made use of the tensile region of the dependence 
of the critical stress rf with the nonglide stress a (cf. eq. (21)). Of course, this dependence is essential to 
characterize the tension/compression asymmetry customarily observed in bcc crystals, cf. Section 1. How¬ 
ever, this is only a weak limitation, as the present CP formulation is sufficiently flexible to admit a full 
(nonlinear) fit to the data shown in Fig. 4. Finally, we emphasize that the present study focuses on plastic 
yielding, and consequently, we have not explored the evolution of the flow stress much beyond the extent 
needed to define a robust yield criterion (cf. Section 3.1). However, this does not detract from the validity of 
the dislocation density evolution model presented in Section 2.4, which has been used prolifically in many 
CP studies (cf. Section 1), and which is being investigated in ongoing studies. 


Acknowledgments 

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore 
National Laboratory under Contract No. DE-AC52-07NA27344. J. M. acknowledges support from DOE’s 
Early Career Research Program. D. C. acknowledges support from the Consejo Social and the PhD program 
of the Universidad Politecnica de Madrid. The authors are indebted to Dr. M. Victoria for inspiration, 
encouragement, and guidance throughout this work. 


^Particularly on non-deviatoric components. 


19 





T[K] 

(b) [110] loading 


Figure 7: Yield strength of W single crystals at the conditions used by Argon & Maloof (1966) in tensile deformation tests under two 
different loading orientations. The experimental data is shown for comparison. The inset shows the results of CP calculations with 
different contributions of the projection tensor activated. 
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Figure 8: Yield strength of W single crystals under the conditions used by Brunner (2010) in uniaxial tensile tests. The experimental 
data is shown for comparison. The inset shows the results of CP calculations with different contributions of the projection tensor 
activated. 
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Figure 9: Contour maps of the yield strength from uniaxial tensile test simulations for 231 uniformly distributed crystallographic 
orientations in the standard triangle at different temperatures. Note that each map has its own distinct numerical scale to aid in the 
visualization of hard and soft regions. 
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Figure 10: Temperature path of the softest and hardest yield directions on the standard triangle as a function of strain rate. 



(a) 



Figure 11: (a) Stress-strain relations at three different strain rates and T = 300 K for a [001] loading orientation, (b) Temperature 
dependence of the yield strength for the softest and hardest directions as a function of strain rate. 
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Figure 12: Dependence of yield strength with strain rate for loading along direction [101] as a function of temperature. The inset 
represents the dependence of the strain rate sensitivity exponent m with temperature. 
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Figure 13: Yield curve for loading along directions ly = [111] and =[11 2] as a function of temperature. 
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Figure 14: Yield surface at 300 K for biaxial loading along directions belonging to the [111] zone. By symmetry, only the 60°-arc need 
be explored. 
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Appendix A. slip systems and latent hardening matrix considered for bcc W. 


Table A.3: Slip systems considered in our calculations, listing the non-normalized crystallographic vectors m^,n^ and n^. Note that in 
DAMASK each slip system is taken both in its positive and negative sense, which is equivalent to formulations where 24 positive slip 
systems are employed (Groger et al., 2008b). 


a 

Reference system 


n" 


1 

[111] (Oil) 

[111] 

[Oil] 

[101] 

2 

[111] (Oil) 

[111] 

[Oil] 

[110] 

3 

[111] (Oil) 

[111] 

[Oil] 

[110] 

4 

[111] (Oil) 

[111] 

[Oil] 

[101] 

5 

[111] (101) 

[111] 

[101] 

[110] 

6 

[111] (101) 

[111] 

[101] 

[Oil] 

7 

[111] (101) 

[111] 

[101] 

[Oil] 

8 

[111] (101) 

[111] 

[101] 

[110] 

9 

[111] (110) 

[111] 

[110] 

[Oil] 

10 

[111] (110) 

[111] 

[110] 

[101] 

11 

[111] (110) 

[111] 

[110] 

[101] 

12 

[111] (110) 

[111] 

[110] 

[Oil] 


Table A.4: Interaction coefficients for the 12 slip systems defined in Table A.3. The letter coding employed is A’: self; ‘CP’: coplanar; 
‘CL’: collinear; ‘O’: orthogonal; ‘G’: glissile; ‘S’: sessile. The reader is referred to Table 1 for the numerical value of each coefficient. 
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Appendix B. Details on the atomistic calculations of non-Schmid parameters 

Critical stresses are computed by applying shear stresses incrementally to a simulation box containing a 
screw dislocation lying on a glide plane forming an angle x with the MRSS plane. The system is schemat¬ 
ically shown in Fig. 3. The box dimensions vary slightly with orientation, such that, for x = the box 
contains 3024 atoms and the dimensions are 21a x 24b x Ic, where a, b, and c are the moduli of the bcc 
lattice vectors x = [121], y = [101], and z = [111], respectively. The calculations are performed using the 
nudged elastic band (NEB) method (Henkelman et al., 2000) implemented in the parallel molecular dy¬ 
namics code LAMMPS (Plimpton, 1995). Periodic boundary conditions are applied along the dislocation 
line direction z while non-periodic and shrink-wrapped boundary conditions are applied along the y and 
X directions. The transition path selected for the NEB calculations is a linear trajectory along the reaction 
coordinate joining two consecutive Peierls valleys, where the dislocation is relaxed to equilibrium. 

Three different forces are applied to different groups of atoms in the simulation box in order to calculate 

These forces recreate the stress tensor (20) in the simulation box: 

1. First, an external force is added to the atoms on the boundary surfaces of the simulation box per¬ 
pendicular to the y-axis to study the T-AT asymmetry. The external force per atom is where 

T is the desired shear stress, is the number of atoms in each nonperiodic surface along z and 

is the cross-sectional area of the each of the bounding surfaces along to y. 

2. To study the contribution from nonglide stresses, an external force is added to the atoms on the 
boundaries of the simulation box perpendicular to the x-axis. The external force per atom is obtained 

fx = where a is the applied nonglide stress, is the number of atoms in each surface and 

LyLz is the cross-sectional area of the each of the surfaces along x. 

3. Further, an external force fy is added to the atoms on the surfaces along the y direction, additionally 
to the shear stress t. fy is defined as fy = with Ny = and is the area of the each of the 
surfaces perpendicular to z. 

31 intermediate replicas are used in the NEB calculations to capture the trajectory and measure the critical 
stress. 
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Appendix C. Crystal plasticity calculations of flow stress dependence with orientation and tempera¬ 
ture 

To demonstrate the performance of the model in the flow stress regime, we carry out calculations for a 
few selected orientations and temperatures up to 10% strain. Figure C.15 shows the stress-strain response at 
a strain rate of 10“^ s“^ as a function of temperature for the [001] loading orientation. This is an orientation 
conducive to multi-slip and thus the system is expected to harden in accordance with eqs. (9), (10), and 
(28) as the deformation progresses. The figure shows results for the full non-Schmid model. For general 
viscoplastic materials it is common to represent the a-e relation as a power law of the type: 

CJ=K£^ (C.l) 

where iC is a constant and n is the so-called hardening exponent. Accordingly, the hardening rate can be 
expressed as: 

^ = Kne^-^ (C.2) 

de 

Fits of eq. (C.l) to the data in Fig. C.15 yield values of n = 0.82, 0.86, and 0.87 for T = 200, 400, and 600 
K, respectively. From eq. (C.2), these numbers result in hardening rates of ^ % 20 MPa/% in all cases. 




Figure C.15: Stress-strain curves for uniaxial loading along the [001] orientation at a strain rate of 10 ^ s ^ for three different 
temperatures. These curves are representative of multi-slip conditions where Taylor-t 3 rpe hardening is enabled. 

Next, we compare the model against the experimental results of Argon & Maloof (1966) for [111] and 
[110] loading at e = 10“^ s“^. To avoid comparing in conditions where twinning may be operative (< 200 
K), which is not captured by our model, we carry out simulations at 293 K. The results are shown in Figure 
C.16, which reveals a good agreement between the full non-Schmid model and the experimental data in 
the [111] loading case. According to Argon & Maloof (1966), yielding under [110] loading occurs at 
approximately 460 MPa, which is immediately followed by an abrupt hardening stage that plateaus at e ^ 
0.8 % to a value of % 760 MPa. Whether or not this is the case, this initial hardening period is not captured 
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by our model. Under both loading conditions, however, the model is seen to reproduce the hardening rates 
in close agreement with the experimental data. 

We emphasize that the results shown in Fig. C.16 have been obtained without fitting to experimental 
(or otherwise) stress-strain curves of any kind, and so the model appears to capture the essential features 
of plastic flow for the selected conditions showcased here. As noted earlier, these preliminary results do not 
imply that the model is suitable for calculating the flow stress under general loading conditions. 
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(a) [111] loading 



(b) [110] loading 

Figure C.16: Flow stress of W single crystals at the conditions used by Argon & Maloof (1966) (cf. Section 3.2) in tensile deformation 
tests under two different loading orientations. The experimental data is shown for comparison. The inset shows the results of CP 
calculations with different contributions of the projection tensor activated. 
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